AD-A191  286 


BJOUtfljUUUMUtaftJUWUiaflWWfllUUUUI  JOBOBOangBBCBGlWnroHinnnwiiroi«nB?»gTCWWl»H»gifmfUgvr8v»M 


AFATL-TR-87-57 


£  &*<=>  /  C=>  oS" 

OTIC  FILE  COE i 


Two-  and  Three-Dimensional 
Computational  Analyses  of 
Projectile-Concrete  Impact 


G  R  Johnson 
R  A  Stryk 


HONEYWEL,  INC 

ARMAMENT  SYSTEMS  DIVISION 
7225  NORTHLAND  DRIVE 
BROOKLYN  PARK,  MINNESOTA  55428 


DECEMBER  1987 


DTIC 


ELECTE 
JAN  2  0  1988 


n 


H 


FINAL  REPORT  FOR  PERIOD  SEPTEMBER  1983-JULY  1987 
I 


APPROVED  FOR  PUBLIC  RELEASE;  DISTRIBUTION  UNLIMITED 


AIR  FORCE  ARMAMENT  LABORATORY 

Air  Force  Systems  Command  I  United  States  Air  Force!  Eglin  Air  Force  Base,  Florida 

88  1  12  162 


, 

l 

! 


t 

> 

k 

j 

q 


NOTICE 


When  Government  drawings,  specifications,  or  other  data  are  used  for 
any  purpose  other  than  in  connection  with  a  definitely  Government-related 
procurement,  the  United  States  Government  incurs  no  responsibility  nor  any 
obligation  whatsoever.  The  fact  that  the  Government  may  have  formulated 
or  in  any  way  supplied  the  said  drawings,  specifications,  or  other  data, 
is  not  to  be  regarded  by  implication  or  otherwise  in  any  manner  construed, 
as  licensing  the  holder,  or  any  other  person  or  corporation;  or  as  conveying 
any  rights  or  permission  to  manufacture,  use,  or  sell  any  patented  invention 
that  may  in  any  way  be  related  thereto. 

The  Public  Affairs  Office  has  reviewed  this  report,  and  it  is  releas¬ 
able  to  the  National  Technical  Information  Service  (NTIS),  where  it  will  be 
available  to  the  general  public,  including  foreign  nationals. 

This  technical  report  has  been  reviewed  and  is  approved  for  publication. 


I 

50 


FOR  THE  COMMANDER 


JCHN  A.\  PALVER,  Col  ,  LEAF 
ChieF,  Munitions  Division 


Even  though  th’s  report  may  contain  special  release  rights  held  by  the 
controlling  office,  please  do  not  request  copies  from  the  Air  Force  Armament 
Laboratory.  If  you  qualify  as  a  recipient,  release  approval  will  be  obtained 
from  the  originating  activity  by  DTIC.  Address  your  request  for  additional 
copies  to: 

Defense  Technical  Information  Center 
Cameron  Station 
Alexandria,  VA  22304-6145 

If  your  address  has  changed,  if  you  wish  to  be  removed  from  our  mailing 
list,  or  if  your  organization  no  longer  employs  the  addressee,  please  notify 
AFATL/M1W,  Eglin  AFB  FL  32542-5434,  to  help  us  maintain  a  current  mailing 
list. 

Copies  of  this  report  should  not  be  returned  unless  return  is  required 
by  security  considerations,  contractual  obligations,  or  notice  on  a  specific 
document. 


nrvwuwjwcwoir.-: 


UNCLASSIFIED  ,  A  „ _ _ 

IfeduAitV  ttAssiFicATi6h  dt  fni*  FaST  /nj>  -  ///7/'z-g^ 

REPORT  DOCUMENTATION  PAGE 

77  report  security  classification  ib.  restrictive  markings 

Unclassified 

it  SECURITY  CLASSIFICATION  AUTHORITY  3.  DISTRIBUTION/ A  VAILABl1 


form  Approved 
OMB  No.  0704-0188 


it  report  security  classification 
Unclassified 


2b  DECLASSIFICATION /DOWNGRADING  SCHEDULE 

4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

7-87-15 

It.  NAME  OF  PERFORMING  ORGANIZATION 

Honeywell 

6b.  OFFICE  SYMBOL 
(If  tppllctble) 

6c  AOORESS  (City,  Sttte.  tnd  ZIP  Code) 

Honeywell  Inc.,  Armament  Systems  Division 

7225  Northland  Drive 

8*  NAME  OF  FUNDING /SPONSORING 
ORGANIZATION 

Munitions  Division 

8b.  OFFICE  SYMBOL 
(If  tppllctble) 

AFATL/MNW 

8c  AOORESS  (City.  Sttte.  tnd  ZIP  Code) 

Air  Force  Armament  Laboratory 
Eglin  AFB,  FL  32542-5434 

6 

3.  DISTRIBUTION /AVAILABILITY  OF  REPORT 

Approved  for  public  release; 
distribution  is  unlimited. 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 

AFATL-TR-87-57 


7«.  NAME  OF  MONITORING  ORGANIZATION 

Clusters  and  Warheads  Branch 
Munitions  Division 


7b.  ADDRESS  (C/ty,  Stitt,  tnd  ZIP  Code) 

Air  Force  Armament  Laboratory 
Eglin  AFB  ,  FL  32542-5434 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 

F 08635- 83-C- 0506 


10.  SOURCE  OF  FUNDING  NUMBERS 


PROJECT 

TASK 

NO. 

NO 

2502 

07 

17 

COSATI  COOES  | 

FIELO 

GROUP 

SUB-GROUP  1 

1 1 .  TITLE  (Include  Security  Cltssrficttion) 

Two-  and  Three-Dimensional  Computational  Analysis  of  Projectile-Concrete  Impact 


1Z  PERSONAL  AUTHOR(S) 

G.R.  Johnson,  R.A.  Stryk 


13a  TYPE  OF  REPORT  13b.  TIME  COVERED  14.  DATE  OF  REPORT  (Yetr,  Month,  Dty)  15.  PAGE  COUNT 

Final  fromScp  83  to  Jul  87  December  1987  55 


16  SUPPLEMENTARY  NOTATION 

Availability  of  this  report  is  specified  on  verso  of  front  cover. 


18.  SUBJECT  TERMS  (Conf/nua  on  reverse  if  necessity  tnd  identify  by  6 lock  number) 

Concrete  Impact,  EPIC 
Hydrocodes  Numerical  Analysis 


19  ABSTRACT  (Confinv*  on  reverse  if  necesstry  tnd  identify  by  block  number) 

This  report  examines  various  computational  approaches  for  projectile  impact  onto  concrete 
targets.  Included  are  two-  and  three-dimensional  computations  performed  with  both  the 
EPIC-2  and  EPIC-3  codes.  Various  techniques,  such  as  tunneling,  eroding,  and  NABOR 
nodes,  are  used  to  compare  results  and  demonstrate  capabilities.  Results  are  presented 
for  normal,  oblique,  and  yawed  impact. 

This  report  also  documents  the  NABOR  algorithms  for  axisymmetric  and  three-dimensional 
geometry.  These  NABOR  algorithms  allow  for  variable  nodal  connectivity  and  are  well  suited 
for  projectile-concrete  impact  problems. 


JO  DISTRIBUTION  /AVAILABILITY  OF  ABSTRACT  21.  ABSTRACT  SECURITY  CLASSIFICATION 

□  unclassifieoiunlimited  iff  same  as  rpt  □  otic  users  Unclassified 


JJ«  NAME  OF  RESPONSIBLE  INDIVIDUAL  J2b.  TELEPHONE  (Include  Are*  Code)  22c  OFFICE  SYMBOL 

MICHAEL  E.  NIXON _ _  (904)  882-2145  AFATL/MNW 


DO  Form  1473,  JUN  86  Previous  editions  ere  obsolete.  SECURITY  CLASSIFICATION  OF  this  page 

Unclassified 


PREFACE 


This  report  was  prepared  by  Honeywell  Inc.,  Armament  Systems  Division,  7225 
Northland  Drive,  Brooklyn  Park,  Minnesota  55428,  for  the  U.S.  Air  Force  Armament 
Laboratory  (AFATL),  Eglin  Air  Force  Base,  Florida  32542,  Under  Contract 
F08635-83-C-0506. 

This  effort  was  conducted  from  September  1983  to  July  1987.  The  authors  would  like  to 
thank  Lts.  Paul  L.  Thee  and  Dennis  L.  May,  AFATL/MNW  program  managers,  and 
William  H.  Cook  and  Michael  E.  Nixon,  AFATL/MNW,  for  many  helpful  technical 
discussions.  Jack  G.  Dodd,  professor  at  Colgate  University  and  consultant  to  Honeywell, 
also  contributed  to  this  work. 

This  report  documents  only  the  portion  of  the  contract  work  concerned  with 
projectile-concrete  impact.  This  contract  was  also  responsible  for  the  development 
of  the  1986  version  of  the  EPIC-2  Code  (Reference  1)  and  the  1987  version  of  the 
EPIC-3  Code  (Reference  2). 


!  L>.  t.:  •jt.ur/ 

;  *v»»  1  >  sbJ  i  tty 

.IVili 

01  s'.  1  i-'-  '•  >- 


iii/iv  (Blank) 


0 1  s  i. 


I 


TABLE  OF  CONTENTS 


Section  Title 

I  INTRODUCTION . 

n  COMPUTATIONAL  STUDY . 

1.  Two-Dimensional  Computations . 

2.  Three-Dimensional  Computations . 

III  NABOR  ALGORITHMS . 

1.  Axisymmetric  Geometry . 

a.  Strains  and  Strain  Rates . 

b.  Deviator  and  Shear  Stresses . 

c.  Pressure  and  Artificial  Viscosity . 

d.  Concentrated  Forces . 

e.  Equations  of  Motion . 

f.  Searching  Algorithm . 

2.  Three-Dimensional  Geometry . 

IV  CONCLUSIONS  AND  RECOMMENDATIONS 

REFERENCES . 


v/vi  (Blank) 


LIST  OF  FIGURES 


Figure 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19 


Title 


Page 


Definition  of  Projectile  Geometry . 3 

EPIC-2  Computation  with  Tunneling  in  Target . 4 

EPIC-2  Computation  with  Erosion  in  Thrget .  5 

EPIC-2  Computation  with  NABOR  Nodes  in  Target .  7 

EPIC-2  Computation  with  Coarse  Grid  and  Erosion  in  Target . 8 

Time  Histories  of  Projectile  Penetration  Depth  and  Velocity  for  the 

EPIC-2  Computations .  9 

Time  Histories  of  Equivalent  Stress  in  the  Projectile  Case  for  the 
EPIC-2  Computations . 

Time  Histories  of  Pressure  in  the  Projectile  Explosive  for  the 
EPIC-2  Computations . 

Time  Histories  of  Axial  Acceleration  at  the  Aft  End  of  the  Projectile 


Time  Histories  of  Projectile  Penetration  Depth  and  Velocity  for  EPIC-2 
and  EPIC-3  Computations . 

Time  Histories  of  Equivalent  Stress  in  the  Projectile  Case  for  the  EPIC-3, 

Normal  Impact  Computations . 

Time  Histories  of  Pressure  in  the  Projectile  Explosive  for  the  EPIC-3, 

Normal  Impact  Computations . 

Time  Histories  of  Axial  and  Transverse  Acceleration  at  the  Aft  End 
of  the  Projectile  for  the  EPIC-3,  Normal  Impact  Computations . 

EPIC-3  Computation  for  Oblique  Impact  (30  Degrees  from  Normal) 
with  Erosion  in  Target . 

EPIC-3  Computation  for  Yawed  Impact  (5  Degrees)  with  Erosion  in  Target ....  24 

Time  Histories  of  Equivalent  Stress  in  the  Projectile  Case  for  the  EPIC-3, 

Oblique  and  Yawed  Impact  Computations .  25 


10 

% 

•iv 

11 

13 

14 

15 

m 

17 

'  iip 

18 

HsilHSEi 

19 

a 

20 

m 

22 

i 

23 

i 

vu 


^%WlV.VkVm\Vm^VA,lV.WA-.Vl'Vl 


LIST  OF  FIGURES  (CONCLUDED) 


Figure  Title  Page 

20  Time  Histories  of  Pressure  in  the  Projectile  Explosive  for  the  EPIC-3, 

Oblique  and  Yawed  Impact  Computations . 26 

21  Schematic  Representation  of  the  NABOR  Computational  Technique  in 

Plane  Strain  Geometry . 28 

22  Definition  of  Velocities  and  Strain  Rates  for  Axisymmetric  Geometry . 29 

23  Volumetric  Strains  for  One-Dimensional  Compression . 32 

24  Plastic  Flow  Relationships  for  Various  Strain  Rates . 33 

25  Definition  of  Nodal  Forces  for  Axisymmetric  Geometry . 36 

26  Derivation  of  Force-Stress  Relationship  for  Hexagonal  Nodal  Arrangement 

in  Axisymmetric  Geometry . 38 

27  Schematic  Representation  of  the  Searching  Algorithm . 39 

28  Arrangement  of  Three-Dimensional  NABOR  Nodes . 42 

29  Attachment  of  Three-Dimensional  NABOR  Nodes  to  Traditional 

Tetrahedral  Element  Grid . 43 


1 


SECTION  I 

INTRODUCTION 


When  a  steel  projectile  impacts  a  concrete  target,  some  special  computational  capabilities 
are  required.  Generally,  the  projectile  will  undergo  only  mild  deformations  and  its 
response  is  best  represented  by  a  Lagrangian  code.  The  concrete  target  often  experiences 
severe  deformations,  which  can  present  problems  for  a  Lagrangian  code. 

It  is  possible  to  predetermine  a  sliding  interface  at  the  axis  of  symmetry  for  normal 
impact  of  a  pointed  projectile.  This  approach  was  used  by  Thigpen  (Reference  3),  who 
used  the  Lagrangian  TOODY  code  to  compute  the  response  of  a  projectile  into 
rock  targets. 

If  the  projectile  is  not  pointed,  then  it  may  not  be  possible  to  accurately  predetermine  the 
sliding  interfaces  in  the  same  manner.  An  alternate  approach,  used  by  Osborn  and 
Matuska  (Reference  4),  is  to  perform  the  computation  in  two  phases.  The  first  phase 
uses  a  rigid  projectile  and  an  Eulerian  target,  and  the  second  phase  uses  a  Lagrangian 
projectile  with  stresses  applied  at  the  boundary.  These  applied  stresses  are  taken  from 
the  first  phase  of  the  computation  and  vary  with  position  and  time.  The  computations 
by  Osborn  and  Matuska  were  performed  with  the  Eulerian  HULL  code  and  the 
Lagrangian  TOODY  code. 

Only  a  limited  number  of  three-dimensional  computations  have  been  performed  to  date. 
Kimsey,  Jonas,  Zukas  and  Johnson  (Reference  5)  used  the  EPIC-3  code  to  determine 
ricochet  conditions  for  a  projectile  impacting  a  concrete  target.  In  these  computations, 
the  concrete  target  material  was  allowed  to  erode,  such  that  it  was  not  necessary  to 
predetermine  the  sliding  interfaces.  This  technique  will  be  described  in  more 
detail  later. 

More  recently,  Rosinsky  (Reference  6)  used  the  DYNA3D  code  to  perform 
three-dimensional,  oblique  impact  computations  for  a  pointed  projectile.  For  these 
computations,  it  was  necessary  to  predetermine  the  penetration  path  by  inserting  a  small 
hole  in  the  target.  Because  the  path  cannot  generally  be  accurately  predicted  at  the 
beginning  of  the  computation,  it  is  usually  necessary  to  perform  an  iterated  series  of 
computations  when  using  this  approach. 

This  report  presents  two-  and  three-dimensional  computations  performed  with  the  EPIC-2 
and  EPIC-3  codes.  An  eroding  target  approach  is  considered,  as  well  as  a  NABOR  node 
approach  that  uses  variable  connectivity.  It  begins  by  evaluating  the  approaches  in 
two-dimensional,  axisymmetric  geometry,  and  then  demonstrates  them  in  three 
dimensions.  Descriptions  of  the  axisymmetric  and  three-dimensional  NABOR  algorithms 
are  also  included. 
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SECTION  II 

COMPUTATIONAL  STUDY 


This  section  presents  a  series  of  two-  and  three-dimensional  computations  performed  with 
the  EPIC-2  and  EPIC-3  codes  (References  1  and  2).  Figure  1  shows  details  of  the 
projectile  used  for  these  computations. 

The  two-dimensional  EPIC-2  grid  uses  triangular  elements  in  a  crossed  triangle  (four 
triangles  per  quadrilateral)  arrangement.  The  three-dimensional  EPIC-3  grid  is 
composed  of  tetrahedral  elements  and  is  similar  to  the  two-dimensional  grid.  In  both 
instances,  sliding  is  allowed  to  occur  between  the  4340  steel  case  and  the  inert  explosive. 
It  is  assumed  that  there  is  no  friction  at  the  steel-explosive  interface  or  at  the 
steel-concrete  interface. 

The  strength  model  for  the  4340  steel  was  developed  by  Johnson  and  Cook  (Reference  7), 
and  the  concrete  model  was  developed  by  Matuska,  Durrett  and  Osborn  (Reference  8). 

For  the  computations  in  this  paper,  the  density  of  the  concrete  is  2308  kg/m3,  and  the 
unconfined  compressive  strength  is  34  megapascal  (MPa).  The  inert  explosive  has  a 
density  of  1817  kg/m3  and  a  yield  stress  of  14  MPa. 

1.  TWO-DIMENSIONAL  COMPUTATIONS 

Figure  2  shows  an  axisymmetric,  tunneling  computation,  where  the  centerline  nodes 
in  the  concrete  target  are  designated  as  slave  nodes,  and  the  outer  surface  of  the  steel 
projectile  is  designated  as  the  master  surface.  Due  to  the  axisymmetric  nature  of  this 
problem,  it  is  possible  to  predetermine  the  path  of  the  projectile,  and  this  allows  for  the 
proper  predetermination  of  the  sliding  interfaces.  From  a  computational  viewpoint,  this 
problem  is  well  defined,  and  it  will  be  the  baseline  to  which  the  other  computations  will 
be  compared. 

Figure  3  shows  an  axisymmetric,  eroding  computation,  where  all  of  the  nodes  in  the 
center  portion  of  the  target  are  designated  as  slave  nodes.  Furthermore,  the  elements  in 
the  target  are  allowed  to  erode  (or  totally  fail)  when  the  equivalent  plastic  strain  exceeds 
1.5.  When  an  element  is  eroded,  it  cannot  develop  any  stresses  or  pressures.  It 
essentially  disappears,  except  that  the  mass  is  retained  at  the  nodes. 

Even  though  this  erosion  approximation  may  appear  to  be  very  severe,  the  depth  of 
penetration  for  this  computation  is  only  11  percent  greater  than  for  the  tunneling 
(baseline)  computation.  Additional  comparisons  will  be  made  later. 

A  distinct  advantage  of  the  erosion  computation  is  that  it  requires  significantly  less 
Central  Processor  Unit  (CPU)  time.  The  erosion  computation  in  Figure  3  required 
only  18  percent  of  the  CPU  time  used  for  the  tunneling  computation  of  Figure  2.  The 
reason  for  this  large  difference  is  that  the  center  line  target  elements  in  the  tunneling 
computation  become  very  small  in  cross  section  as  they  move  radially  outward  along  the 
outside  of  the  case.  The  small  cross  section  leads  to  a  small  minimum  altitude,  which  in 
turn  leads  to  a  small  integration  time  increment.  The  net  result  is  that  the  tunneling 
computation  requires  many  more  integration  cycles. 
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The  reason  the  time  increment  remains  much  larger  for  the  erosion  computation  is  that 
the  center  line  elements  are  eroded  when  they  become  highly  strained,  and  the  eroded 
elements  cannot  govern  the  integration  time  increment.  This  leaves  lower-strained 
elements,  with  larger  altitudes,  to  govern  the  time  increment. 

Another  distinct  advantage  of  the  eroding  target  approach  is  that  it  is  not  necessary  to 
predetermine  the  path  of  the  projectile  or  to  predetermine  the  precise  sliding  interfaces. 
This  is  especially  important  in  three  dimensions.  All  of  the  nodes  in  the  appropriate 
portion  of  the  target  are  designated  as  slave  nodes,  and  the  projectile  simply  follows  the 
path  of  the  least  resistance,  as  governed  by  the  mechanics  in  the  numerical  algorithm. 

Figure  4  shows  an  axisymmetric  NABOR  computation,  where  the  nodes  in  the  center 
portion  of  the  target  are  designated  as  NABOR  nodes.  They  are  also  designated  as  slave 
nodes.  The  NABOR  option  allows  for  variable  nodal  connectivity.  Each  node  is  affected 
only  by  its  nearest  neighbor  nodes.  As  the  nodes  move  closer  than  their  equilibrium 
distance,  they  generate  compressive,  repulsive  forces.  Conversely,  when  they  move  apart 
they  can  generate  tensile  attractive  forces.  Material  strength  effects  are  also  included. 
The  key  to  this  approach  is  that  it  is  possible  to  have  variable  nodal  connectivity.  A 
node  can  take  on  new  nearest  neighbors  thus  allowing  all  forms  of  distortion.  Based  on 
the  concept  of  nearest  neighbors,  this  option  has  been  designated  the  NABOR  option. 

The  specific  algorithm  will  be  presented  in  the  next  section. 

The  NABOR  results  shown  in  Figure  4  are  in  good  general  agreement  with  the  previous 
computations,  with  the  depth  of  penetration  being  14  percent  greater  than  the  depth 
obtained  from  the  tunneling  (baseline)  computation. 

Some  other  features  of  this  approach  should  also  be  noted.  Because  it  is  a  Lagrangian 
algorithm  (the  mass  moves  with  the  grid),  it  is  possible  to  connect  the  NABOR  nodes  to 
the  traditional  finite  element  grid  as  shown  in  Figure  4.  It  is  also  possible  to  designate 
the  NABOR  nodes  as  slave  nodes  such  that  they  can  slide  along  the  master  surface  on 
the  projectile.  Like  the  eroding  target  approach  discussed  previously,  there  is  no  need  to 
predetermine  the  penetration  path  or  the  precise  sliding  interface.  The  path  is 
automatically  determined  from  the  mechanics  in  the  numerical  algorithm. 

Although  the  results  of  the  eroding  computation  in  Figure  3  appear  to  be  similar  to  the 
results  of  the  NABOR  computation  in  Figure  4,  the  two  methods  may  not  be  similar  for 
all  cases.  The  eroding  target  approach  may  not  be  as  accurate  as  the  NABOR  approach 
for  blunt  nose  shapes.  Conversely,  the  NABOR  approach  incorporates  many  simplifying 
assumptions  and  its  accuracy  has  not  been  verified  for  a  wide  range  of  problems. 
Nevertheless,  the  initial  results  with  both  approaches  are  very  encouraging. 

Another  factor  that  should  be  considered  when  performing  comparative  computations  is 
the  effect  of  the  grid  size.  Figure  5  shows  the  results  of  a  coarser  grid  in  the  target 
when  using  the  eroding  target  approach.  Here,  the  number  of  target  elements  is  only 
half  as  many  as  used  in  the  computations  of  Figures  2  and  3.  When  compared  to  the 
eroding  target  computation  of  Figure  3,  the  coarser  grid  computation  of  Figure  5 
produces  a  maximum  penetration  depth  which  is  only  3  percent  less. 

Figure  6  shows  penetration  depth  and  velocity  for  the  four  computations  shown  in  Figures 
2  through  5.  The  tunneling  (baseline)  computation  is  probably  the  most  accurate  but  the 
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others  give  similar  results.  The  NABOR  computation  provides  the  most  deceleration 
during  the  initial  0.003  second,  but  it  then  decelerates  more  slowly  than  the  others,  and 
eventually  experiences  the  greatest  depth  of  penetration. 

Another  comparison  is  provided  in  Figure  7,  where  the  equivalent  stress  at  a  selected 
element  in  the  steel  case  is  shown  as  a  function  of  time.  Generally,  the  responses  are 
very  similar  and  all  four  cases  show  a  small  amount  of  equivalent  plastic  strain.  The 
two  eroding  computations  show  more  oscillation  than  the  tunneling  computation,  and 
this  is  probably  due  to  the  instantaneous  changes  in  the  resistance  of  the  elements  that 
erode.  The  NABOR  results  show  a  higher  stress  at  later  times,  and  this  occurs  because 
the  projectile  is  still  decelerating  as  indicated  in  Figure  6. 

A  similar  comparison  is  shown  in  Figure  8  where  the  pressure  in  the  inert  explosive  is 
shown  as  a  function  of  time.  Again,  the  general  responses  are  very  similar,  with  the 
tunneling  (baseline)  computation  having  less  oscillation  than  the  other  responses. 

The  final  comparison  is  shown  in  Figure  9,  where  the  acceleration  at  the  top  of  the 
projectile  case  is  shown  as  a  function  of  time.  Due  to  the  high  frequency  responses,  it  is 
difficult  to  make  meaningful  comparisons.  The  maximum  rigid  body  deceleration  is  only 
3000-4000  Gs.  For  the  tunneling  and  the  two  eroding  computations,  there  appears  to  be 
a  lower  frequency  in  the  range  of  0.005  to  0.010  second.  During  this  intermediate  time 
range,  the  inert  explosive  is  compressed  and  loses  contact  with  the  top  of  the  case.  This 
allows  the  top  of  the  case  to  vibrate  as  a  plate,  and  this  probably  leads  to  the  lower 
frequency  response.  The  subsequent  higher  frequency  response  is  probably  due  to  the 
inert  explosive  rebounding  back  against  the  end  plate  of  the  case. 

The  NABOR  computation  provides  a  noisier  signal  for  all  of  the  time  history  responses 
in  Figures  7  through  9.  The  reason  for  this  trend  is  not  obvious. 


2.  THREE-DIMENSIONAL  COMPUTATIONS 

Although  the  two-dimensional  axisymmetric  computations  provide  for  meaningful 
comparisons  between  various  computational  approaches,  the  three-dimensional,  oblique, 
and  yawed  impact  problems  are  of  the  most  interest.  The  three-dimensional  model  for 
the  projectile  is  shown  in  Figure  1  and  the  complete  model  for  the  erosion  computations 
is  shown  in  Figure  10,  The  three-dimensional  model  for  a  target  with  NABOR  nodes  is 
similar,  with  the  exception  that  it  is  rectangular  in  shape. 

Figure  11  shows  a  three-dimensional  computation  of  the  normal  impact  problem  that  was 
previously  analyzed  in  axisymmetric  geometry.  All  of  the  nodes  in  the  center  portion  of 
the  target  are  designated  as  slave  nodes,  and  the  target  elements  are  allowed  to  erode  at 
an  equivalent  plastic  strain  of  1.5.  There  is  no  predetermination  of  penetration  path  or 
sliding  interfaces. 

As  shown  in  the  figure,  the  projectile  does  not  travel  in  a  straight  line,  as  expected,  but 
instead  rotates  6  degrees.  The  initial  off-axis  perturbation  is  due  to  the  sliding  interface 
algorithm.  If  a  3lave  node  in  the  target  is  exactly  positioned  on  a  line  connecting  two 
triangular  master  surfaces,  the  equations  of  motion  are  adjusted  in  a  direction  normal  to 
only  one  of  the  master  surfaces.  In  the  initial  geometry,  there  are  numerous  occurrences 


Figure  9.  Time  Histories  of  Axial  Acceleration  at  the  Aft  End  of  the  Projectile  for 
the  EPIC'2  Computations 
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of  this  match  up  between  slave  nodes  and  master  surface  interfaces.  Although  the  initial 
disturbance  can  be  explained,  the  magnitude  of  the  eventual  rotation  is  surprising. 

The  maximum  penetration  depth  result  is  in  excellent  agreement  with  the  two-dimensional, 
axisymmetric  results.  It  is  only  9  percent  greater  than  the  tunneling  (baseline)  depth  of 
Figure  2,  and  it  falls  between  the  two  depths  obtained  with  the  eroding  target 
computations  of  Figures  3  and  5. 

Figure  12  shows  a  three-dimensional  NABOR  computation  of  the  same  normal  impact 
problem.  Here,  there  is  a  small  projectile  rotation  of  3  degrees,  but  this  rotation  is  much 
less  than  that  experienced  in  the  erosion  computation  of  Figure  11.  The  penetration 
depth  is  similar  to  that  of  the  other  computations,  being  only  1  percent  less  than  the 
tunneling  (baseline)  depth  of  Figure  2. 

Penetration  depths  and  velocities  for  the  three-dimensional  computations  of  Figures  11 
and  12  are  shown  in  Figure  13.  As  was  the  case  for  the  two-dimensional  NABOR 
computations,  the  three-dimensional  NABOR  computations  also  decelerate  the  projectile 
at  a  slightly  higher  rate  than  do  the  other  computations.  Generally,  the 
three-dimensional  computations  (eroding,  NABOR)  show  good  agreement  with  the 
two-dimensional  tunneling  (baseline)  computation. 

Figure  14  shows  equivalent  stresses  in  selected  elements  as  a  function  of  time.  These 
elements  correspond  to  the  element  monitored  in  the  two-dimensional  computations,  as 
shown  in  Figure  7.  The  responses  at  points  A  and  D  should  be  identical  if  the  responses 
are  truly  symmetrical.  For  the  eroding  computation,  there  is  a  noticeable  difference 
between  the  responses  at  the  two  locations,  and  this  is  consistent  with  the  overall 
projectile  rotation  as  shown  in  Figure  11.  There  appears  to  be  more  consistency  between 
the  two  locations  for  the  NABOR  computation,  and  it  is  probably  because  the  projectile 
does  not  experience  as  much  rotation  as  it  does  in  the  eroding  computation. 

Also,  a  small  amount  of  plastic  strain  is  experienced  at  all  locations  and  it  is  slightly 
higher  for  the  NABOR  computations.  This  is  most  likely  caused  by  the  higher 
deceleration  as  shown  in  Figure  13. 

Although  there  are  some  definite  differences  between  the  two-dimensional  results  of 
Figure  7  and  the  three-dimensional  results  of  Figure  14,  the  equivalent  stress  responses 
are  generally  very  comparable. 

Figure  15  shows  pressure  as  a  function  of  time  for  the  same  location  monitored  in  the 
two-dimensional  computations.  Here,  the  responses  are  very  similar  at  points  B  and  E, 
and  they  do  not  appear  to  be  noticeably  affected  by  the  rotation  of  the  projectile.  The 
pressures  for  the  three-dimensional  NABOR  computation  are  slightly  higher  than  for  the 
corresponding  eroding  computation,  and  this  is  probably  due  to  the  higher  deceleration. 
Again,  however,  the  two-dimensional  results  of  Figure  8  and  the  three-dimensional 
results  of  Figure  15  are  generally  comparable. 
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Figure  16  shows  acceleration  as  a  function  of  time  at  the  aft  end  of  the  projectile.  When 
the  three-dimensional  results  of  Figure  16  are  compared  to  the  two-dimensional  results  of 
Figure  9,  there  appears  to  be  a  distinct  difference.  The  magnitude  of  the  three-dimensional 
accelerations  is  much  lower  than  that  experienced  in  the  two-dimensional  computations. 

A  possible  explanation  is  that  the  three-dimensional  grid  is  coarser  (it  has  no  center 
nodes  as  does  the  two-dimensional  crossed  triangles)  and  cannot  develop  as  high  a 
frequency  response  as  can  the  two-dimensional  grid.  Also,  a  transverse  acceleration  is 
developed  from  the  asymmetric  response  of  the  projectile. 

The  three-dimensional,  normal  impact  computations  of  Figures  11  and  12  have  been 
performed  to  evaluate  the  adequacy  of  the  three-dimensional  computing  techniques.  The 
primary  problems  of  interest  involve  oblique  and/or  yawed  impact,  such  as  shown  in 
Figures  17  and  18.  Both  of  these  computations  use  the  eroding  target  approach,  but  it 
would  be  expected  that  the  NABOR  approach  would  give  similar  results. 

The  oblique  impact  computation  shows  the  projectile  experiencing  a  significant  rotation 
of  19  degrees,  from  an  initial  orientation  of  30  degrees  to  a  final  orientation  of  49  degrees. 
Similarly,  the  yawed  impact  computation  goes  from  an  initial  orientation  of  5  degrees  to 
a  final  orientation  of  16  degrees,  for  a  total  rotation  of  11  degrees. 

The  corresponding  equivalent  stresses  and  pressures  are  shown  as  a  function  of  time  in 
Figures  19  and  20.  The  equivalent  stress  at  point  D  rises  faster  than  at  point  A.  This 
occurs  because  the  compressive  axial  stress  and  the  compressive  bending  stress  combine 
at  this  point  during  early  time.  Although  the  details  of  these  responses  vary  from  those 
of  normal  impacts,  in  some  respects  they  are  very  similar.  The  equivalent  plastic  strains, 
for  instance,  are  not  significantly  different  from  those  experienced  during  normal  impact. 

The  pressures  are  very  similar  to  those  experienced  during  normal  impact,  with  the 
pressures  at  point  B  being  slightly  higher  than  those  at  point  E. 
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Figure  20.  Time  Histories  of  Pressure  in  the  Projectile  Explosive  for  the  EPIC-3, 
Oblique  and  Yawed  Impact  Computations 
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This  section  presents  the  specific  algorithms  for  both  the  axisymmetric  and  the 
three-dimensional  NABOR  techniques.  These  are  the  algorithms  used  for  the 
computations  in  Figures  4  and  12  in  the  preceding  section.  These  are  Lagrangian 
algorithms  because  the  material  moves  with  the  grid.  The  distinguishing  feature  of 
these  approaches  is  that  they  allow  for  variable  nodal  connectivity. 

The  NABOR  technique  was  first  demonstrated  in  plane  strain  geometry  during  the  early 
phases  of  this  contract  (Reference  9).  Based  on  those  results,  it  was  then  extended  to 
axisymmetric  and  three-dimensional  geometry.  It  should  be  noted  that  the  idea  of 
variable  nodal  connectivity  is  not  new.  A  Particle-and-Force  (PAF)  method  was 
developed  for  fluids  in  the  early  1960’s  by  Daly,  Harlow  and  Welch  (Reference 
10).  Although  there  are  some  similarities  in  the  general  approaches  of  the  PAF  and 
NABOR  algorithms,  they  are  quite  different. 

A  schematic  representation  of  a  plane  strain  NABOR  grid  for  an  impact  problem  is 
shown  in  Figure  21.  The  NABOR  nodes  are  somewhat  analogous  to  circular  flexible 
disks.  The  equilibrium  position  is  a  hexagonal  close-pack,  where  each  internal  node  has 
six  neighbors.  As  node  j  moves  closer  to  node  i,  a  compressive,  repulsive  force  is  exerted 
on  node  i,  as  shown  in  the  lower  portion  of  Figure  21.  (This  force  also  acts  on  node  j.) 
Node  k  is  moving  away  from  node  i  and,  therefore,  exerts  a  tensile,  attractive  force  on 
node  i.  Variable  connectivity  is  also  demonstrated  in  Figure  21,  where  node  m  is  moving 
inward  to  become  a  neighbor  of  node  i,  and  node  j  is  moving  outward  to  eventually  have 
no  effect  on  node  i. 

1.  AXISYMMETRIC  GEOMETRY 

The  axisymmetric  geometry  is  similar  to  the  plane  stress  geometry,  with  the 
complicating  exceptions  that  the  hoop  strains,  strain  rates,  and  stresses  must  all  be 
considered.  The  mass  of  a  NABOR  node  for  axisymmetric  geometry  is 

M  =  n  V3  r^Do2  (1) 

where  r0  is  the  initial  radial  coordinate  at  the  center  of  the  node,  p0  is  the  density  of 
the  material,  and  D0  is  the  initial  diameter  of  all  NABOR  nodes  contained  in  the 
specific  problem. 


a.  Strains  and  Strain  Rates 

Figure  22  shows  a  bond  between  NABOR  nodes  i  and  j.  D  is  the  current 
distance  between  nodes,  and  Uj  and  vt  are  the  r  and  z  velocities  of  node  i.  V,N  is 
the  normal  velocity  vector  along  the  bond  at  node  i,  and  V,p  is  the  velocity  vector 
perpendicular  to  the  bond.  .The  following  equations  represent  various  strains  and 
strain  rates. 
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Figure  2 1 .  Schematic  Representation  of  the  NABOR  Computational  Technique  in 
Plane  Strain  Geometry 
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Figure  22.  Definition  of  Velocities  and  Strain  Rates  for  Axisymmetric  Geometry 
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The  strain  rates  along  the  bond,  perpendicular  to  the  bond,  and  in  the  hoop  direction 
are  given  by  £  N  and  £  p  and  £  e,  respectively.  While  £  N  and  £  9  can  be  exactly  defined, 

£  p  must  be  approximated. 

It  will  later  be  shown  that  ft  is  often  within  the  range,  --1  <  p  <  0,  and  it  is  therefore 
assumed  that  fi  =  -0.5.  The  shear  strain  rate,  y  ,  is  also  an  approximation  inasmuch  as 
three  nodes  (rather  than  two)  are  required  to  give  an  exact  description. 

The  equivalent  strain  rate,  £ ,  can  be  integrated  to  give  an  approximate  equivalent 
strain,  £.  In  Equation  7,  At  is  the  integration  time  increment,  and  a  is  a  factor  that 
varies  linearly  from  0  at  D  =  1.3D0  to  1.0  for  D  <  1.15D'0.  This  a  factor  is  also  applied 
to  stresses  and  pressures;  it  provides  for  a  smooth  transition  at  D  =  1.3D0  and  lessens 
the  influence  of  nodes  at  extended  bond  distances. 

N  is  the  number  of  bonds  (neighbor  nodes)  at  the  NABOR  node,  where  a  bond  is  defined 
as  two  nodes  with  a  distance  between  them  of  D  <  1.3D'„.  The  other  term  is 
D'0  =  D0  vr0/r.  It  is  the  effective,  unstrained  nodal  diameter  in  the  r-z  plane,  caused  by 
the  current  radial  coordinate,  r,  being  different  from  the  initial  radial  coordinate,  r0. 

Although  £  is  updated  from  a  bond  strain  rate,  £  ,  it  must  be  carried  from  cycle  to  cycle 
as  a  nodal  quantity.  The  increment  of  nodal  strain  is,  therefore,  the  average  of  all  the 
incremental  bond  strains  associated  with  that  node.  Because  the  nodal  connectivity  can 
vary  from  cycle  to  cycle,  there  can  be  no  retention  of  bond  variables;  only  the  nodes  can 
retain  variables  from  cycle  to  cycle. 

The  volumetric  strain  for  pressure  computations  is  m  =  p/p0  -  1  =  V0/V  -  1,  where 
p0  and  p  are  the  initial  and  current  densities  of  a  specified  element  of  material,  and  V0 
and  V  are  the  corresponding  volumes.  Because  a  two-dimensional  vomme  cannot  be 
determined  from  the  distance  between  two  nodes,  p  is  determined  from: 

Mi  =  D'0/D  -  1.1547  (8a) 

Mz  =  (D0/D)2  -  1  (8b) 


M  =  max  (mi,  M2) 


(8c) 
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Equation  8a  insures  that  two  nodes  cannot  get  too  close  together.  The  constant  in  this 
equation  (1.1547)  allows  pure  shear  to  occur  without  introducing  a  significant  volumetric 
strain.  This  also  reduces  the  locking  and  cracking  problems  associated  with  the  initial 
plain  strain  NABOR  algorithm  (Reference  9).  In  equation  8b,  D  =  ED  (D0/D'0)/N  is  the 
average  factored  distance  to  all  of  the  neighbor  nodes  with  D  <  1.3D'0.  Similarily,  D'0  is 
the  average  of  D'0  for  the  two  nodes  forming  the  bond. 

Figure  23  shows  the  necessity  for  p2.  The  initial  hexagonal  geometry  is  uniformly 
strained  in  the  vertical  direction  only.  Furthermore,  it  is  assumed  that  D'0  =  D0  for  all 
nodes.  Therefore,  \i  should  be  identical  for  all  bonds.  This  is  clearly  not  the  case  for  p,; 
the  distance  between  nodes  i  and  j  is  greater  than  the  distance  between  nodes  i  and  k. 
The  computed  p,  is  incorrect  for  both  cases.  Fortunately,  the  nodal  quantity,  n  *  cr/h,  is 
essentially  correct  for  small  volumetric  strains.  In  equation  8c,  p2  is  the  average  of  the 
two  nodal  volumetric  strains. 

Another  geometric  distance  that  will  be  used  later  for  the  force  computations  is 

Dp  =  (D'0)2/D  (1  +  M)  (9) 

where  Dp  is  a  representative  nodal  distance  in  the  direction  perpendicular  to  the  bond 
distance,  D. 

b.  Deviator  and  Shear  Stresses 

The  deviator  and  shear  stresses  are  dependent  on  the  strains  and  strain  rates  in 
Equations  2-7.  Because  the  strains  and  strain  rates  are  based  on  approximations,  the 
deviator  and  shear  stresses  are  also  approximately  determined.  Only  plastic  flow  is 
considered  and  elastic  effects  are  not  included.  The  deviator  stresses  in  the  normal 
(along  the  bond),  and  hoop  directions,  and  the  shear  stress  are  given  by 


sN  =  2aae^/3£ 

(10) 

s9  =  2 a  a  e0/3£ 

(11) 

r  =  a  5  y  /3i 

(12) 

where  a  is  the  equivalent  tensile  flow  stress,  and  eN  and  e,  are  deviator  strain 
rates.  The  equivalent  tensile  flow  stress  in  the  EPIC  codes  (References  1  and  2)  is 
expressed  as  a  function  of  the  strain,  strain  rate,  temperature  and/or  pressure. 

The  deviator  and  shear  stresses  are  dependent  on  strain  rates,  which  have  been 
approximated  in  some  instances.  Assuming  i9  =  o,  the  effect  of  setting  p  =  -0.5  in 
Equation  3  is  shown  in  Figure  24.  For  incompressible  plastic  flow,  £P  =  -tN  so  that 
P  =  -1.0.  On  the  other  extreme,  p  ~  0  for  one-dimensional  wave  propagation  along 
tN.  Fortunately,  sN  is  not  very  sensitive  to  p  in  this  range,  and  the  assumption  that 
P  =  -0.5  is  acceptable. 
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c.  Pressure  and  Artificial  Viscosity 


The  total  normal  stress  along  the  bond  and  the  hoop  stress  are 

<tn  =  sN  -  (P  +  Q)  (13) 

ae  =  s9  -  (P  +  Q)  (14) 

where  sN  and  se  are  the  previously  defined  deviator  stresses,  P  is  the  hydrostatic  pressure 
and  Q  is  the  artificial  viscosity.  Compression  is  positive  for  P  and  Q.  The  artificial 
viscosity  damps  out  localized  oscillations  that  would  otherwise  occur  during  wave 
propagation  (Reference  11). 


Q  =  CL  p  csD|£ v|  +  CqpD2£V2 
Q  =  0 


for  i\  <  0  (15) 

for  £V  >  0 


where  CL  and  Cq  are  dimensionless  coefficients,  c,  is  the  sound  velocity,  £v  is  the 
volumetric  strain  rate  and  the  other  terms  are  as  previously  defined. 

The  pressure  is  determined  from  the  Mie-Gruneisen  equation  of  state  (Reference  12).  For 
compressive  strains  (p  >  0) 

P  =  (K,p  +  K2p2  +  K3p3)  (1  -  Tp/2)  +  I>E  (16) 

where  Kj,  K2,  K3  and  T  are  material-dependent  constants  and  E  is  internal  energy  per 
unit  mass.  At  high  pressures,  the  pressure  can  be  significantly  affected  by  the  internal 
energy.  Therefore,  it  is  necessary  to  solve  the  pressure  and  energy  equations 
simultaneously.  Equation  16  can  be  rewritten  as 

F  +  At  =  PH  +  rpE‘  +  At  (17) 

where 

PH  =  (Kjp  +  K2p2  +  K3p3)  (1  -  I>/2)  (18) 

E‘  +  At  =  El  +  AE,  +  (FADN  +  F,rAR,N,  +  FfARjNJ)/(M,  +  Mj)  (19) 

The  superscripts  t  and  t  +  At  represent  the  times  at  the  beginning  and  end  of  an 
integration  cycle.  In  Equation  19,  AE3  is  the  increment  of  internal  energy  per  unit  mass 
generated  by  the  shear  and  deviator  stresses  during  the  current  integration  cycle.  In  the 
same  equation,  N  is  the  average  number  of  neighbor  nodes  of  the  two  nodes  forming  the 
bond,  F  is  the  average  force  (from  P  and  Q  only)  acting  at  that  bond,  and  AD  =  D£NAt  is 
the  distance  through  which  the  force  acts.  The  other  forces  (F,R,  Ff),  distances  (AR„  AR,), 
and  numbers  of  neighbor  nodes  (Njf  Nj)  are  for  the  hoop  components.  The  masses  of  the 
two  nodes  are  M,  and  Mj. 


The  average  forces  (F,  F?,  Ff )  are  functions  of  (P  +  Q)‘  and  (P  +  Q)1  *  4t.  The  specific 
force-stress  relationships  will  be  presented  later.  Because  (P  +  Q)‘  is  known  from  the 
previous  cycle,  and  Q‘  +  At  is  known  from  the  current  cycle,  the  only  unknown  on  the 
right  side  of  Equation  19  is  P  +  At.  Substituting  Equation  19  into  Equation  17  gives  an 
explicit  equation  for  Pl  +  4t  that  is  consistent  with  E‘  +  At. 

There  are  two  significant  assumptions  associated  with  N  and  F.  They  are  related  to  the 
fact  that  a  single  node  usually  has  more  than  one  neighbor,  and  therefore,  more  than 
one  bond.  When  this  occurs,  it  is  assumed  that  the  other  bonds  generate  the  same 
amount  of  energy  as  the  specific  bond  being  considered.  This  assumption  is  performed  by 
the  presence  of  N. 

The  next  assumption  concerns  F.  In  the  force-stress  relationships,  (P  +  Q)1  is  a  nodal 
quantity,  which  is  the  average  of  (P  +  Q)  for  all  the  individual  bonds  at  time  =  t. 

As  mentioned  previously,  it  is  not  possible  to  carry  bond  variables  from  cycle  to  cycle 
because  the  nodal  connectivity  is  not  fixed.  Therefore,  the  nodal  quantity  (P  +  Q)1  will 
not  be  the  exact  bond  quantity  at  time  =  t.  This  approximation  tends  to  reduce 
(P  +  Q)1  at  the  bonds  where  the  pressure  is  the  highest  and  where  the  most  internal 
energy  is  generated.  As  a  result  of  this  reducing  assumption,  the  computed  internal 
energy  is  sometimes  too  low,  and  this  leads  to  a  computed  energy  loss  (internal  and 
kinetic)  for  the  entire  system. 

For  expanded  states  (p  <  0),  K2  =  K3  =  0  and  the  maximum  negative  pressure  is 
limited  to  a  fraction  (about  half)  of  the  strength  of  the  material,  a.  Furthermore,  the 
tensile  pressure  (P  +  Q)  is  multiplied  by  a.  This  gives  a  smooth  transition  of 
a(P  +  Q)  =  0  at  D  =  1.3D0. 

d.  Concentrated  Forces 

The  concentrated  forces  acting  on  the  nodes  are  statically  equivalent  to  the  normal, 
hoop,  and  shear  stresses  (on,  as,  t).  The  comparable  normal,  hoop,  and  shear  forces 
(Fn,  F9,  Fp)  are  shown  in  Figure  25  and  are  expressed  as 


Fn  =  2nf  Dpc7N/V3" 

(20) 

Fe  =  27t<7*A/N 

(21) 

Fp  =  2rrrDpT/v/3_ 

(22) 

In  Equations  20  and  22,  r  is  the  average  radial  coordinate  of  the  two  nodes  forming 
the  bond,  and  Dp  is  the  perpendicular  distance  given  in  Equation  9.  In  Equation  21, 

A  =  v,3(D'0)2/2  is  the  cross-sectional  area  of  the  NABOR  node,  and  N  is  the  number  of 
bonds  on  the  node.  This  allows  a  portion  of  the  hoop  force  to  be  computed  for  every 
bond.  This  hoop  calculation  is  performed  separately  for  each  of  the  two  nodes. 
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These  relationships  are  derived  for  a  hexagonal  close-pack  geometry  as  shown  in 
Figure  26.  Unfortunately,  for  other  geometric  arrangements,  the  force-stress  relationship 
is  different.  For  example,  a  checkerboard  arrangement  of  rows  and  columns,  where 
there  are  only  four  neighbors,  gives  FN  =  DpaN.  Although  the  work  performed 
to  date  has  used  only  the  relationships  of  Equations  20-22,  it  may  eventually  be 
preferable  to  alter  these  relationships  such  that  they  become  a  function  of  the  current 
geometric  arrangement. 

After  the  normal,  hoop,  and  shear  forces  (FN,  Fg,  Fp)  have  been  determined,  the  r  and  z 
components  (Fr,  F,)  are  determined  as  shown  in  Figure  25.  For  node  i 

Fr‘  =  — Fjf  cos0  -  Fp  sin0  -  Fg  (23) 

Fi  =  -Fn  sin<A  +  Fp  cos<A  (24) 

The  forces  on  node  j  (with  the  exception  of  Fg)  are  equal  and  opposite  to  those  on  node  i. 
For  each  specific  bond,  there  is  a  static  moment  unbalance  of  D*Fp.  These  localized 
unbalances  tend  to  cancel  one  another  and  appear  to  have  little  affect  on  the  overall 
system  response. 

e.  Equations  of  Motion 

The  nodal  equations  of  motion  for  node  i,  the  r  direction,  are 

ur  =  uf"  +  EF’At/Mi  (25) 

r-  +  41  =  rf  +  ut+At  (26) 

where  u‘+  is  the  constant  velocity  for  the  time  increment  between  t  and  t  +  At,  and  u‘~ 
is  the  velocity  prior  to  time  =  t.  The  net  force  in  the  r  direction  is  EF,!,  which  contains 
contributions  from  all  bonds  containing  node  i.  It  can  also  accept  contributions  from 
traditional  triangular  and/or  quad  elements,  such  as  demonstrated  in  Figure  4.  The 
average  integration  time  increment  (about  time  =  t)  is  At,  and  Mi  is  the  mass  of  node  i. 

Equation  26  gives  the  updated  nodal  positions  r-  *  4t.  The  integration  time  increment, 

At,  must  be  less  than  the  minimum  sound  speed  transit  time  between  any  two  nodes. 

The  equations  of  motion  in  the  z  direction  have  the  same  form  as  Equations  25  and  26. 

f.  Searching  Algorithm 

A  schematic  of  the  searching  algorithm  is  shown  in  Figure  27.  The  neighbor  nodes 
of  NABOR  node  i  are  all  the  nodes  within  a  distance  D  <  1.3  D0  to  node  i.  The 
neighborhood  nodes  all  have  distances  D  >  1.3  D0,  but  they  are  the  next  closest  nodes  to 
NABOR  node  i.  The  total  number  of  neighbor  nodes  plus  neighborhood  nodes,  in  this 
example,  is  eight.  For  an  interior  NABOR  node  in  a  hexagonal  close  pack  geometry, 
there  would  be  six  neighbor  nodes  and  two  neighborhood  nodes. 
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Figure  26.  Derivation  of  Force-Stress  Relationship  for  Hexagonal  Nodal  Arrange¬ 
ment  in  Axisymmetric  Geometry 
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Only  the  neighbor  nodes  can  affect  node  i  directly  by  transmitting  forces  to  it.  It  is 
helpful  to  monitor  the  neighborhood  nodes  because  they  are  the  nodes  most  likely  to 
become  neighbor  nodes.  The  important  features  of  the  searching  algorithm  are  as  follow. 

For  every  NABOR  node,  the  distances  to  all  neighbor  and  neighborhood  nodes  are 
computed  for  each  cycle  of  integration.  This  allows  neighborhood  nodes  to  become 
neighbor  nodes  as  soon  as  D  <  1.3  D'0. 

Every  NABOR  node  has  a  timer  variable,  T„  which  is  the  minimum  time  required  for  a 
new  node  (other  than  a  neighbor  or  neighborhood  node)  to  possibly  become  a  new 
neighbor.  It  is  expressed  as 

T,  =  (Dy  -  D'0)/Vref  (27) 

where  D,j  is  the  distance  between  two  nodes  i  and  j,  and  V^f  is  a  reference  velocity 
given  by 

Vref  =  max(V0,  AVy)  (28) 

Here,  V0  is  a  user-supplied  input  velocity,  and  AVtj  is  the  computed  velocity  difference 
(magnitude  only)  between  nodes  i  and  j. 

If  V0  is  larger  than  any  relative  closing  velocity  experienced  during  the  course  of  the 
event,  then  every  neighbor  node  would  be  found  before  it  gets  too  close  to  the  associated 
NABOR  node.  This  may  lead  to  excessive  searching,  especially  for  relatively  high 
velocity  impact  problems. 

At  the  other  extreme,  if  V0  is  very  small,  the  Vref  will  be  governed  by  the  computed 
velocity  differences,  AVy.  This  would  reduce  the  searching,  and  every  new  neighbor  node 
would  probably  be  found  before  it  gets  too  close  to  the  associated  NABOR  node.  There 
may  be  some  instances  when  the  new  neighbor  nodes  would  not  be  found  soon  enough. 
When  this  occurs,  it  is  necessaiy  to  restart  the  computation  from  an  earlier  time,  using  a 
higher  value  of  V0.  For  every  cycle  of  integration,  T,  is  reduced  by  At. 

Tt  *  “  =  T\  -  At  (29) 

When  Ts  becomes  negative,  a  global  search  is  initiated.  This  search  examines  every 
other  NABOR  node  in  the  grid  and  establishes  an  updated  set  of  neighbor  and 
neighborhood  nodes.  It  also  establishes  an  updated  T,  based  on  distances  and  relative 
velocities  of  all  nodes,  except  the  neighbor  and  neighborhood  nodes. 

It  is  at  this  point  that  the  value  of  monitoring  the  neighborhood  nodes  can  be  seen.  If 
the  neighborhood  nodes  were  not  monitored,  the  T,  would  often  be  governed  by  one  of  the 
neighborhood  nodes,  and  it  would  often  provide  a  much  smaller  time  until  the  next 
search,  thus  causing  more  searching. 
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For  a  limited  number  of  computations  performed  to  date,  the  computing  time  required 
for  searching  has  not  been  excessive.  It  is  evident  from  Equations  27  and  28  that  higher 
impact  velocities  cause  more  searching.  A  counteracting  feature  is  that  the  time 
increment  does  not  significantly  decrease  for  the  NABOR  computations  because  the 
minimum  distance  between  nodes  is  always  D  =  D'0.  For  high  impact  velocities,  using 
traditional  Lagrangian  techniques,  the  integration  time  increment  often  becomes  very 
small  due  to  the  highly  distorted  grid.  This  results  in  more  integration  cycles  (and 
increased  computing  time)  for  the  traditional  Lagrangian  approach. 

2.  THREE-DIMENSIONAL  GEOMETRY 

The  three-dimensional  NABOR  algorithm  follows  the  same  general  approach  as  the 
preceding  axisymmetric  algorithm.  In  some  aspects  it  is  more  complicated,  due  to  the 
third  dimension;  but  in  other  respects  it  is  less  complicated,  because  all  of  the 
hoop-related  complications  of  axisymmetric  geometry  are  eliminated. 

This  subsection  includes  a  description  of  the  three-dimensional  mass  distribution, 
the  strains  and  strain  rates,  and  the  concentrated  forces.  Other  portions  of  the 
three-dimensional  algorithm  are  very  similar  to  those  of  the  axisymmetric  algorithm 
and  will  not  be  repeated.  This  includes  the  deviator  and  shear  stresses,  the  pressure 
and  artificial  viscosity,  the  equations  of  motion,  and  the  searching  algorithm. 

The  initial  NABOR  node  geometry  consists  of  nested  layers  of  nodes,  where  each  layer 
has  a  hexagonal  close-pack  arrangement,  as  shown  in  Figure  28.  The  projected 
hexagonal  area  is  V3D02/2,  and  the  effective  layer  thickness  is  V6D0/3.  The  resulting 
NABOR  node  mass  is 

M  =  (30) 

Figure  29  shows  how  the  NABOR  nodes  are  attached  to  the  traditional  tetrahedral 
element  grid.  Here,  the  layers  of  nodes  are  in  the  x-z  plane.  This  basic  grid 
arrangement  was  used  for  the  three-dimensional  NABOR  computation  in  Figure  12. 

The  three  dimensional  volumetric  strains  are  computed  in  a  manner  similar  to  that  used 
in  two  dimensions. 


Mi  =  D0/D  -  1.2247 

(31a) 

M2  =  (D0/D)3  -  1 

(31b) 

M  *  max  (pi,  m2) 

(31c) 

where  Dn  is  the  initial  NABOR  node  diameter,  D  is  the  current  distance  between  nodes, 
and  D  is  the  average  distance  of  all  neighbor  nodes  that  have  bonds  (D  <  1.3D0)  with 
the  NABOR  node.  The  constant  in  equation  31a  (1.2247)  allows  pure  shear  to  occur 
without  introducing  a  significant  volumetric  strain. 


Figure  29.  Attachment  of  Three-Dimensional  NABOR  Nodes  to  Traditional  Tetra¬ 
hedral  Element  Grid 


The  normal  and  shear  strain  rates  in  three  dimensions  are  similar  to  the 
two-dimensional  strain  rates  in  Figure  22  and  in  Equations  2  and  5.  The  specific 
formulation  for  the  normal  strain  rate  between  nodes  i  and  j  is 


/N  =  (A*AVX  +  B*AVy  +  C*AVZ)/D  (32) 

where  A,  B,  and  C  are  direction  cosines  from  node  i  to  node  j,  and  the  velocity  differences 
(AV,,  AVy,  AVj)  have  the  form  AVX  =  Uj  -  Ui.  The  velocities  in  the  x,  y,  and  z  directions 
are  designated  as  u,  v,  and  w,  respectively. 

There  are  two  additional  normal  strain  rates  that  are  perpendicular  to  the  bond. 

Because  these  cannot  be  determined  directly,  they  are  estimated  to  be 

£  Pi  =  £  P2  =  -0.25 £N  (33) 

This  is  essentially  the  same  as  setting  p  =  -0.5  in  Equation  4. 

The  net  three-dimensional  shear  strain  rate  is  the  net  relative  velocity  (between  nodes 
i  and  j)  that  lies  in  a  plane  perpendicular  to  the  bond  (line  connecting  the  two  nodes), 
divided  by  the  distance  between  the  two  nodes.  It  is  expressed  as 


y  =  vy  *  +  Yy  +  y 


where 


Y  x  =  [(B2  +  C2)AVx  -  A*B*AVy  -  A.C.AVJ/D  (35) 

Y  y  =  [-A*B«AVX  +  (A2  +  C2)AVy  -  B.C-AVJ/D  (36) 

y\  =  (-A-C-AV,  -  B*C*AVy  +  (A2  +  B2)AVJ/D  (37) 

The  equivalent  strain  rate,  i  ,  and  equivalent  strain,  £,  are  as  presented  previously  in 
Equations  6  and  7. 

After  the  strains  and  strain  rates  are  defined,  the  stresses  and  pressures  can  be 
determined  in  a  manner  similar  to  that  used  in  Equations  10  - 19.  The  resulting  stresses 
of  interest  are  the  total  normal  stress  along  the  bond,  <jn,  and  the  shear  stress,  t. 

The  lower  right  portion  of  Figure  28  shows  how  the  stresses  can  be  converted  into 
forces.  The  resulting  normal  and  shear  forces  are  expressed  as 

Fn  =  V2Dp2<tn/4  (38) 

F,  =  V2D2x/4 


where  Dp2  =  D03/D(l  +  p)  tends  to  adjust  the  projected  cross-sectional  area  normal  to 
the  bond,  based  on  the  volumetric  strain  and  the  distance  between  the  nodes. 


Now,  the  x,  y,  and  z  forces  on  node  i  can  be  obtained. 

Fx‘  =  A-Fn  +  (yx/y)Fs 
Fy  =  B*Fn  +  (y  y/y  )F, 

F?  =  C«Fn  +  (y  Jy  )FS 

The  forces  on  node  j  axe  equal  and  opposite  to  those  on  node  i. 


The  equations  of  motion  and  the  searching  algorithm  are  similar  to  those  of  th6 
two-dimensional  case. 


SECTION  IV 

CONCLUSIONS  AND  RECOMMENDATIONS 


This  report  has  presented  a  series  of  two-  and  three-dimensional  computations  for  a  steel 
projectile  impacting  a  concrete  target.  In  two  dimensions,  it  was  demonstrated  that  the 
eroding  target  and  NABOR  target  approaches  could  be  used  to  give  results  that  agreed 
very  well  with  the  more  exact  tunneling  (baseline)  approach.  In  three  dimensions,  the 
normal  impact  conditions  were  repeated,  and  although  there  were  some  discrepancies, 
the  results  are  in  general  agreement  with  the  two-dimensional  results.  Examples  of 
oblique  and  yawed  impacts  were  also  presented  to  demonstrate  the  three-dimensional 
capability.  For  these  two  example  problems,  there  were  no  significant  differences  in  the 
monitored  equivalent  stresses  or  pressures  when  compared  to  the  normal  impact  results. 

The  axisymmetric  and  three-dimensional  NABOR  algorithms  have  also  been  presented. 
The  initial  results  are  very  encouraging,  but  more  computations  must  yet  be  performed 
to  thoroughly  assess  this  technique. 
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